% Computes the steady state of twoSbank.mod
% File sends "SS values" to twoSbank.mod and updates the parameters
% that represent some steady state ratios (= "steady state parameters")
% Christoph Gortz

function [ys,check] = twosb2v8CQ48_steadystate(ys,exe)
global M_

%% 1. I load the values of the deep parameters in a loop.
%%
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
end              

check = 0; % check=0 limits steady state values to real numbers (see  Dynare user guide page 71).


%% 2. Compute the values for the steady state parameters
%%
% Imposed
Lss = 1;
cbeta = 1/(1+constebeta/100);
ga = cga/100 ;
gv = cgv/100 ;
clandap = 1 + clandapaux;
clandaw = 1 + clandawaux;
crho = 1/(cminrho - 1); % restricting cminrho in (0,1)
ga = cga/100;
spreadss = 0.005; 


% Values of steady state parameters needed for financial intermediaries:
% Steady state of bank's stochastic discount factor
sdfbss = 1;
% Steady state nominal interest rate
rss = 1/cbeta * exp(ga+(alfac)/(1-alfai)*gv) * (constepinfc/100 + 1);
% "private" leverage ratio varrho => q*s/n = lrdata = 1/(1-psi)*varrho
lriss = (1-psiiss)*lridata;
lrcss = (1-psicss)*lrcdata;

% varomega is chosen so that risk premium is 40 basis points per quarter in SS
varomegai = (1-cthetabss*(spreadss*lriss+(rss))*1/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv))*((1-psiiss)*lridata)^(-1);
varomegac = varomegai; % institutional structure of bank should be symmetric

% Resulting value for Rbcss, Rbiss
rbiss = ( ((1-varomegai*(1-psiiss)*lridata)*1/(cthetabss*1/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv)) - rss)*1/lriss + rss ) / grconstepinfc;
rbcss = ( ((1-varomegac*(1-psicss)*lrcdata)*1/(cthetabss*1/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv)) - rss)*1/lrcss + rss ) / grconstepinfc;
% same as: rbciss = (rss+spreadss)/grconstepinfc;



% 2.1 Steady state parameter values for the simple two sector model:

% 2.1.1 Steady state values without adjustment costs a la Huffman and Wynne
% (used for initial guess)

iX1aux = (1-alfai)^(alfai-1)*alfai^-alfai*(1-alfac)^(1-alfai)*alfac^(alfac*(1-alfai)/(1-alfac));
% steady state rental rate
icrk = (( exp((1/(1-alfai))*gv)/cbeta - (1-cdeltai) )*iX1aux*(1/clandap)^((alfac-alfai)/(1-alfac)) )^((1-alfac)/(1-alfai)); 
% steady state relative price of investment
icpi = iX1aux* icrk^((alfai-alfac)/(1-alfac))*(1/clandap)^((alfac-alfai)/(1-alfac));
% steady state wage
icw = ( (1/clandap)*(1-alfac)^(1-alfac)*alfac^alfac*icrk^-alfac )^(1/(1-alfac));
% steady state capital labor ratio in consumption sector
iKcLc = icw/icrk*(alfac/(1-alfac)); 
% steady state capital labor ratio in investment sector
iKiLi = icw/icrk*(alfai/(1-alfai)); 
% Labour
iLc = Lss*( 1+ nvalueic*(iKcLc^alfac)/(iKiLi^alfai)*icpi^-1 )^-1;
iLi = Lss - iLc;
% steady state capital in consumption sector
iKc = iKcLc*iLc;
% steady state capital in investment sector
iKi = iKiLi*iLi;
% steady state investment sector output
iiss = (1/clandap)*iKiLi^alfai*iLi;
% steady state investement in I sector:
iiiss = iKi*(exp(1/(1-alfai)*gv) - (1-cdeltai));
% steady state consumption in I sector:
iicss = iKc*(exp(1/(1-alfai)*gv) - (1-cdeltac));




% 2.2 Calculate values of "steady state parameters"
% Steady state equations for 2 Sector model with adjustment costs a la 
% Huffman and Wynne and financial intermediaries

% Relabeling the variables:
% ki = x(1), kc = x(2), ii = x(3), ic = x(4) rik = x(5) rck = x(6) Li = x(7) Lc = x(8) w = x(9)  pi = x(10)

% % F contains the model equations in SS (all equal to zero):
 F = @(x) [
% % Equation ki
 x(1) - x(9)/x(5) * alfai/(1-alfai) * x(7);
% % Equation kc
 x(2) - x(9)/x(6) * alfac/(1-alfac) * x(8);
% % Equation w 
 x(9) - ( 1/clandap * (1-alfac)^(1-alfac)*alfac^alfac*x(6)^(-alfac) )^(1/(1-alfac));
% % Equation pi
x(10) - ( 1/clandap*(1-alfac)^(1-alfac)*alfac^alfac * ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltac) )^-alfac * ((x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(4)^(-crho-1))^-alfac )  /  (  1/clandap*(1-alfai)^(1-alfai)*alfai^alfai * ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltai) )^-alfai * ((x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(3)^(-crho-1))^-alfai    )^((1-alfac)/(1-alfai));
% % Equation Lc
 x(8) - (1+nvalueic*(x(2)/x(8))^alfac/((x(1)/x(7))^alfai) * x(10)^(-1) )^(-1);
 % % Equation Li
 x(7) - Lss + x(8);
% % Equation rki
 x(5) - ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltai) )*x(10) * (x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(3)^(-crho-1);
 % % Equation rkc
 x(6) - ( (exp(1/(1-alfai)*gv))/cbeta - (1-cdeltac) )*x(10) * (x(3)^(-crho)+x(4)^(-crho))^(-1/crho -1)*x(4)^(-crho-1);
 % % Equation ii
 x(3) - x(1)*(exp(1/(1-alfai)*gv) - (1-cdeltai));
 % % Equation ic
 x(4) - x(2)*(exp(1/(1-alfai)*gv) - (1-cdeltac));
 ];

 x0 = [iKi; iKc; iiiss; iicss; icrk; icrk; iLi; iLc; icw; icpi;];     % Initial guess for start values
 options = optimset('Display','off','TolFun',1e-12,'TolX',1e-12);
 x = fsolve(F,x0,options);  % x contains the variables x(1)-x(10)

X1aux = (1-alfai)^(alfai-1)*alfai^-alfai*(1-alfac)^(1-alfai)*alfac^(alfac*(1-alfai)/(1-alfac));
crki = x(5);
crkc = x(6);
cpi = x(10);
cw = x(9);
KcLc = x(2)/x(8);
KiLi = x(1)/x(7);
Lc = x(8);
Li = x(7);
css = (1/clandap)*KcLc^alfac*Lc;
Kc = x(2);
Ki = x(1);
K = x(2) + x(1);
Kibar = Ki*exp((1/(1-alfai))*gv);
Kcbar = Kc*exp((1/(1-alfai))*gv);
iiss = x(3);
icss = x(4);
iss = (1/clandap)*KiLi^alfai*Li; 




% Solving for the bank's steady state parameters
% Consumption Sector:
% Relabeling the variables:
% nucs = x(1), etacs = x(2), clambc = x(3),
% % F contains the model equations in SS (all equal to zero):
Fx = @(xx) [
% % Equation nu; substituted Equation z2iss combined with equation for leverage ratio
xx(1) - (1-cthetabss)*sdfbss*exp(-ga-alfac/(1-alfai)*gv)*(spreadss) - cthetabss*cbeta* (  ((spreadss)*  xx(2)/(xx(3)-xx(1))   + rss)/grconstepinfc  ) *xx(1);
% Equation eta
xx(2) - (1-cthetabss)*sdfbss*exp(-ga-alfac/(1-alfai)*gv)*rss - cthetabss*cbeta* (  ((spreadss)*  xx(2)/(xx(3)-xx(1))   + rss)/grconstepinfc  ) *xx(2);
 % Leverage ratio for clamb
lrcss - xx(2)/(xx(3)-xx(1));
 ];

% Initial guesses for bank's steady state values:
gnui = 0.5;
getai = 0.5;
gclambi = 1.0;

xx0 = [gnui; getai; gclambi; ];     % Initial guess for start values
options = optimset('Display','off','TolFun',1e-12,'TolX',1e-12);
xx = fsolve(Fx,xx0,options);  % x contains the variables xx(1)-xx(10)
nucss = xx(1); 
etacss = xx(2);
clamb = xx(3);

% For Investment sector:
nuiss = nucss;
etaiss = etacss;

% Steady state of price of capital in I sector
qiss = (iiss^(-crho) + icss^(-crho))^(-1/crho-1) * iiss^(-crho-1) * cpi;
% Steady state of price of capital in C sector
qcss = (iiss^(-crho) + icss^(-crho))^(-1/crho-1) * icss^(-crho-1) * cpi;
% Quantity of financial claims held by bank in both sectors
siss = Kibar;
scss = Kcbar;
% Steady state output
yss = (css+cpi*iss) + ctau*psicss*qcss*scss + ctau*psiiss*qiss*siss;
%yss = (1/(1-cg))*(css+cpi*iss) + ctau*psicss*qcss*scss + ctau*psiiss*qiss*siss;
% Steady state of z2 for Bank in both sectors
z2iss =  ((spreadss)*  lriss + rss)*1/grconstepinfc;
z2css =  ((spreadss)*  lrcss + rss)*1/grconstepinfc;
% Steady state of z1 for Bank in both sectors
z1css = z2css;
z1iss = z2iss;
% Steady state of shock to capital quality in both sectors
sncss = 1;
sniss = 1;





%% 5. Assign SS variables for vector ys
%%
 
labis = 0;
labcs = 0;
labs =0;
kcs = 0;
kis = 0;
kpis = 0;
kpcs = 0;
mccs = 0;
mcis = 0;
zcapis = 0;
zcapcs = 0;
rkis = 0;
rkcs = 0;
pis = 0;
cs = 0;
lams = 0;
phifis = 0;
phifcs = 0;
inveis = 0;
invecs = 0;
inves = 0;
ysss = 0;
pinfcs = 0;
pinfis = 0;
ws =0;
rs = 0;
zs = 0; 
vs = 0; 
z_ls = 0;
v_ls = 0;
bs = 0;
gs = 0;
qss = 0;
mss = 0; 
sws = 0;
ewmas = 0;
%=====================================================================
% NOTE zeros when usinh demeaned data
dys = 0 + 0 + alfac/(1-alfai)*0 + (cga + alfac/(1-alfai)*cgv)*0;  % Observables...
dcs = 0-0 + 0 + alfac/(1-alfai)*0 + (cga + alfac/(1-alfai)*cgv)*0;
dinves = 0-0 + 1/(1-alfai)*0 + ( 1/(1-alfai)*cgv)*0;
dpis = 0 + 0 + (alfac-1)/(1-alfai)*0 + (cga+ (alfac-1)/(1-alfai)*cgv)*0;
dws = 0 + 0 + alfac/(1-alfai)*0 + (cga + alfac/(1-alfai)*cgv)*0;
labobss = 0;
robss = 0 + (log(constepinfc/100 + 1) - log(cbeta) + (cga + alfac/(1-alfai)*cgv))*0;
pinfobscs = 1*(0) + constepinfc*0;
pinfobsis = 1*(0) + constepinfi*0 ;
%pinfobsis = pinfobscs + dpis;
spreadobscs = 0 + (log(rbcss) - robss)*0;
spreadobsis = 0 + (log(rbiss) - robss)*0;
dequitys = exp(ga + alfac/(1-alfai)*gv)*(cga + alfac/(1-alfai)*cgv)*0;
%======================================================================

%spreads = ( 0+log(rbcss)-robss )* ( exp(0+log(qcss))*exp(0+log(scss)) )/exp(0+log(qcss*scss*lrcdata^-1))... 
%        + ( 0+log(rbiss)-robss )* ( exp(0+log(qiss))*exp(0+log(siss)) )/exp(0+log(qiss*siss*lridata^-1));
%dcredits = ( cga+ alfac/(1-alfai)*cgv )* ( exp(0+log(qcss))*exp(0+log(scss)) )/ (  exp(0+log(qcss))*exp(0+log(scss)) + exp(0+log(qiss))*exp(0+log(siss)) )...
%         + (  cga + alfac/(1-alfai)*cgv )* ( exp(0+log(qiss))*exp(0+log(siss)) )/ (  exp(0+log(qcss))*exp(0+log(scss)) + exp(0+log(qiss))*exp(0+log(siss)) );

rrates = 0;     % ... end observables.
epinfmacs = 0;    
epinfmais = 0;
spinfcs = 0;
spinfis = 0;
spis = 0;
spkcs = 0;
spkis = 0;
zcapifs = zcapis;    % flexible economy variables. Should have same ss as above.
zcapcfs = zcapcs;
rkcfs = rkcs;
rkifs = rkis;
kpifs = kpis;
kpcfs = kpcs;
kifs = kis;
kcfs = kcs;
lamfs = lams;
phififs = phifis;
phifcfs = phifcs;
cfs = cs;
invefs = inves;
inveifs = inveis;
invecfs = invecs;
yfs = ysss;
labfs = labs;
labifs = labis;
labcfs = labcs;
pifs = pis;
wfs = ws;
sdfbfs = 0;     % beginning of variables specific for banking sector...
nucfs = 0;      % (flexible price economy)
nuifs = 0;
rbcfs = 0;
rbifs = 0;
z1cfs = 0;
z1ifs = 0;
etaifs = 0;
etacfs = 0;
z2cfs = 0;
z2ifs = 0;
lrcfs = 0;
lrifs = 0;
ncfs = 0;
nifs = 0;
qcfs = 0;
qifs = 0;
scfs = 0;
sifs = 0;   % ...end of variables specific for banking sector
sdfbs = 0;  % beginning of variables specific for banking sector...
rbcs = 0;   % (sticky price economy)
rbis = 0;
z1cs = 0;
z1is = 0;
z2cs = 0;
z2is = 0;
nucs = 0;
nuis = 0;
etacs = 0;
etais = 0;
lrcs = 0;
lris = 0;
qcs = 0;
qis = 0;
scs = 0;
sis = 0;
ncs = 0;
nis = 0;
rpcs = 0;
rpis = 0;   
sncs = 0;
snis = 0;
psics = 0;
psiis = 0; % ...end of variables specific for banking sector
epkc1aux = 0;
epkc2aux = 0;
epkc3aux = 0;
epkc4aux = 0;
epkc8aux = 0;
epki1aux = 0;
epki2aux = 0;
epki3aux = 0;
epki4aux = 0;
epki8aux = 0;
ez1aux = 0;
ez4aux = 0;
ez8aux = 0;
ev1aux = 0;
ev4aux = 0;
ev8aux = 0;
cthetabs = 0;

ys = [ labis; labcs; labs; kcs; kis; kpis; kpcs; mccs; mcis; zcapis; zcapcs; rkis; rkcs;
pis; lams; phifis; phifcs; cs; inveis; invecs; inves; ysss; pinfcs; pinfis; ws; rs; zs;
vs; z_ls; v_ls; bs; gs; qss; mss; sws; ewmas; dys; dcs; dinves;    dpis;     dws; labobss; robss; pinfobscs; 
pinfobsis; spreadobscs; spreadobsis; rrates;      epinfmacs; epinfmais; spinfcs; spinfis; spis; spkcs; spkis; zcapifs; zcapcfs; 
rkcfs; rkifs; kpifs; kpcfs; kifs; kcfs; lamfs; phififs; phifcfs; cfs; invefs; inveifs;
invecfs; yfs; labfs; labifs; labcfs; pifs; wfs; 
sdfbfs; nucfs; nuifs; rbcfs; rbifs; z1cfs; z1ifs; etaifs; etacfs; z2cfs; z2ifs; lrcfs; lrifs; ncfs; nifs; qcfs; qifs; scfs; sifs;
sdfbs; rbcs; rbis; z1cs; z1is; z2cs; z2is; nucs; nuis; etacs; etais; lrcs; lris; qcs; qis; scs; sis; ncs; nis; rpcs; rpis; sncs; snis; psics; psiis;
epkc1aux; epkc2aux; epkc3aux; epkc4aux; epkc8aux; epki1aux; epki2aux; epki3aux; epki4aux; epki8aux;
ez1aux; ez4aux; ez8aux; ev1aux; ev4aux; ev8aux;
cthetabs; dequitys;
]; 

%% 6. update the "steady state parameters" in the global workspace  
%% steady state parameters: X1aux crki crkc cpi cw KcLc KiLi Lc Li css iss
%% Kc Ki K Kibar Kcbar yss iiss icss

id1 = strmatch('X1aux',M_.param_names,'exact');     % Get the index of crki in M_.params
id2 = strmatch('crki',M_.param_names,'exact');      % Get the index of crkc in M_.params
id3 = strmatch('crkc',M_.param_names,'exact');                   
id4 = strmatch('cpi',M_.param_names,'exact');
id5 = strmatch('cw',M_.param_names,'exact');                   
id6 = strmatch('KcLc',M_.param_names,'exact');
id7 = strmatch('KiLi',M_.param_names,'exact');
id8 = strmatch('Lc',M_.param_names,'exact');
id9 = strmatch('Li',M_.param_names,'exact');
id10 = strmatch('css',M_.param_names,'exact');
id11 = strmatch('iss',M_.param_names,'exact');         
id12 = strmatch('Kc',M_.param_names,'exact');           
id13 = strmatch('Ki',M_.param_names,'exact');                   
id14 = strmatch('K',M_.param_names,'exact');
id15 = strmatch('Kibar',M_.param_names,'exact');                   
id16 = strmatch('Kcbar',M_.param_names,'exact');
id17 = strmatch('yss',M_.param_names,'exact');
id18 = strmatch('iiss',M_.param_names,'exact');
id19 = strmatch('icss',M_.param_names,'exact');
id20 = strmatch('z1css',M_.param_names,'exact');
id21 = strmatch('z1iss',M_.param_names,'exact');     
id22 = strmatch('z2css',M_.param_names,'exact');  
id23 = strmatch('z2iss',M_.param_names,'exact');                   
id24 = strmatch('rbcss',M_.param_names,'exact');
id25 = strmatch('rbiss',M_.param_names,'exact');                   
id26 = strmatch('rss',M_.param_names,'exact');
id27 = strmatch('lrcss',M_.param_names,'exact');
id28 = strmatch('lriss',M_.param_names,'exact');
id29 = strmatch('nucss',M_.param_names,'exact');
id30 = strmatch('nuiss',M_.param_names,'exact');
id31 = strmatch('qcss',M_.param_names,'exact');         
id32 = strmatch('qiss',M_.param_names,'exact');           
id33 = strmatch('sncss',M_.param_names,'exact');                   
id34 = strmatch('sniss',M_.param_names,'exact');
id35 = strmatch('scss',M_.param_names,'exact');                   
id36 = strmatch('siss',M_.param_names,'exact');
id37 = strmatch('clamb',M_.param_names,'exact');
M_.params(id1) = X1aux;
M_.params(id2) = crki;
M_.params(id3) = crkc;
M_.params(id4) = cpi;
M_.params(id5) = cw;
M_.params(id6) = KcLc;
M_.params(id7) = KiLi;
M_.params(id8) = Lc;
M_.params(id9) = Li;
M_.params(id10) = css;
M_.params(id11) = iss;
M_.params(id12) = Kc;
M_.params(id13) = Ki;
M_.params(id14) = K;
M_.params(id15) = Kibar;
M_.params(id16) = Kcbar;
M_.params(id17) = yss;
M_.params(id18) = iiss;
M_.params(id19) = icss;
M_.params(id20) = z1css;
M_.params(id21) = z1iss;
M_.params(id22) = z2css;
M_.params(id23) = z2iss;
M_.params(id24) = rbcss;
M_.params(id25) = rbiss;
M_.params(id26) = rss;
M_.params(id27) = lrcss;
M_.params(id28) = lriss;
M_.params(id29) = nucss;
M_.params(id30) = nuiss;
M_.params(id31) = qcss;
M_.params(id32) = qiss;
M_.params(id33) = sncss;
M_.params(id34) = sniss;
M_.params(id35) = scss;
M_.params(id36) = siss;
M_.params(id37) = clamb;


